Systematic study of d-wave superconductivity in the 2D repulsive Hubbard model 
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The cluster size dependence of superconductivity in the conventional two-dimensional Hubbard 
model, commonly believed to describe high-temperature superconductors, is systematically studied 
using the Dynamical Cluster Approximation and Quantum Monte Carlo simulations as cluster solver. 
Due to the non-locality of the d-wave superconducting order parameter, the results on small clusters 
show large size and geometry effects. In large enough clusters, the results are independent of the 
cluster size and display a finite temperature instability to d-wave superconductivity. 
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Despite years of active research, the understanding of 
pairing in the high-temperature "cuprate" superconduc- 
tors (HTSC) remains one of the most important out- 
standing problems in condensed matter physics. While 
conventional superconductors are well described by the 
BCS theory, the pairing mechanism in HTSC is believed 
to be of entirely different nature. Strong electronic cor- 
relations play a crucial role in HTSC, not only for su- 
perconductivity but also for their unusual normal state 
behavior. Hence, models describing itinerant correlated 
electrons, in particular the two-dimensional (2D) Hub- 
bard model and its strong-coupling limit, the 2D t-J 
model, were proposed to capture the essential physics 
of the CuO-planes in HTSC 0,0. Despite the fact that 
these models are among the mostly studied models in 
condensed matter physics, the question of whether they 
contain enough ingredients to describe HTSC remains an 
unsolved problem. 

Many different techniques, from analytic to numerical 
have been applied to study superconductivity in these 
models. The Mermin- Wagner theorem Q and the rig- 
orous results in Ref. Q preclude d x 2_ y 2 superconduct- 
ing long-range order at finite temperatures in the 2D 
models. Superconductivity may however exist - as in 
the attractive Hubbard model - as topological order at 
finite temperatures below the Kosterlitz-Thouless (KT) 
transition temperature |f| . Recent renormalization group 
studies indicate that the ground-state of the doped weak- 
coupling 2D Hubbard model is superconducting with 
a ci c 2_j / 2-wave order parameter The possibility of 
d x 2_ y 2-wave pairing in the 2D Hubbard and t-J models 
was also indicated in a number of numerical studies of fi- 
nite system size (for a review see 

0)- 

Only recent numer- 
ical calculations for the t-J model provided evidence for 
pairing at T = in relatively large systems for physically 
relevant values of J/t ||. Quantum Monte Carlo (QMC) 
simulations are also employed to search for such a tran- 
sition These studies indicate an enhancement of the 
pairing correlations in the d x 2_ y 2 channel with decreasing 
temperature. Unfortunately the Fermion sign problem 
limits these studies to temperatures too high to study a 
possible KT transition. Another difficulty of these meth- 



ods arises from their strong finite size effects, often ruling 
out the reliable extraction of low-energy scales. In fact, a 
reliable finite-size scalin g ha s only recently been achieved 
in the negative-U model [lfj , where the relevant tempera- 
ture scales are much higher. The available results for the 
positive-U model so far have thus been inconclusive, and 
a treatment within a non-perturbative scheme that goes 
beyond the conventional finite size techniques is clearly 
necessary to resolve the controversy as to whether there 
exists finite temperature superconductivity in these mod- 
els. 

In this Letter we use the Dynamical Cluster Approxi- 
mation (DC A) (for a review see 01) to explore the 
superconducting instability in the 2D Hubbard model 
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where cfj (creates) destroys an electron with spin a on 
site i, iiic is the corresponding number operator, t the 
hopping amplitude between nearest neighbors (...) and 
U the on-site Coulomb repulsion. In the DCA we take 
advantage of the short length-scale of spin correlations 
in optimally doped HTSC |l3| to map the original lattice 
model onto a periodic cluster of size N c — L c x L c embed- 
ded in a self-consistent host. Thus, correlations up to a 
range £ < L c are treated accurately, while the physics on 
longer length-scales is described at the mean-field level. 
By increasing the cluster size, it thus allows us to sys- 
tematically interpolate between the single-site dynamical 
mean-field result and the exact result while remaining in 
the thermodynamic limit. We solve the cluster problem 
using QMC simulations H3- 

We present results of large cluster calculations - up to 
26 sites - that indicate that the 2D Hubbard model has a 
superconducting instability at a finite temperature. This 
conclusion is reached due to several factors: Simulations 
on small clusters, where d-wave order is topologically al- 
lowed, show large finite size and geometry effects leading 
to inconclusive results. However, since the average sign 
in DCA QMC simulations is significantly larger than in 
finite-size QMC counterparts, exploring lower tempera- 
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tures and larger clusters becomes possible. In addition, 
the advent of new parallel vector machines, such as the 
CRAY XI at ORNL, improves the speed of these calcula- 
tions by more than one order of magnitude compared to 
conventional architectures, making simulations on large 
clusters with a small average sign feasible. Within the 
limits of current computational capability, we observe 
finite transition temperatures in the largest affordable 
clusters. There the results are independent of cluster 
size within the error bars, although we cannot preclude a 
further small reduction in transition temperatures in yet 
larger clusters. 

Previous DCA simulations with a cluster of four sites, 
the smallest cluster that can capture d x 2_ y 2-wave pair- 
ing, with U equal to the bandwidth W — St, show good 
general agreement with HTSC • In the paramagnetic 
state, the low-energy spin excitations become suppressed 
below the crossover temperature T*, and a pseudogap 
opens in the density of states at the chemical potential. 
At lower temperatures, we find a finite temperature tran- 
sition to antiferromagnetic long-range order at low dop- 
ing, while at larger doping, the system displays an insta- 
bility to d x 2„ J( 2-wave superconducting long-range order. 
This apparent violation of the Mermin- Wagner theorem 
is a consequence of the small cluster size studied (see also 
More recent results obtained with a similar quan- 
tum cluster algorithm confirm the presence of antiferro- 
magnetism and superconductivity in the groundstate of 
the 2D Hubbard model 

With increasing cluster size however, the DCA progres- 
sively includes longer-ranged fluctuations while retaining 
some mean- field character. Larger clusters are thus ex- 
pected to systematically drive the Need temperature to 
zero and hence recover the Mermin- Wagner theorem in 
the infinite cluster size limit. In contrast, superconduc- 
tivity may persist as KT order even for large cluster sizes. 

Since the large cluster simulations presented here are at 
the limit of current computational capabilities, we are re- 
stricted in our ability to explore both the parameter space 
and different cluster sizes. We choose the parameters to 
favor superconducting and antiferromagnetic order. In 
our study of superconductivity, we choose U = At = W/2 
(we take t as our unit of energy). While we observe that 
larger values of U yield higher transition temperatures in 
the 4-site cluster, the smaller value of U greatly reduces 
the sign problem and thus allows us to simulate larger 
cluster sizes. We focus on a doping of 10%, where the 
pairing correlations are maximal for U — W/2. To study 
antifcrromagnetism, we focus on the undoped model and 
set U = 8t, where the Neel temperature is highest. 

Furthermore, we have to be careful in selecting differ- 
ent cluster sizes and geometries. Much can be learned 
from simulations of finite size systems, where periodic 
boundary conditions are typically used. Betts and Flynn 
[lij systematically studied the 2D Heisenberg model on 
finite size clusters and developed a grading scheme to de- 



termine which clusters should be used. The main quali- 
fication is the "imperfection" of the near-neighbor shells: 
a measure of the (in)completeness of each neighbor shell 
compared to the infinite lattice. In finite size scaling cal- 
culations they found that the results for the most perfect 
clusters fall on a scaling curve, while the imperfect clus- 
ters generally produce results off the curve. Here, we 
employ some of the cluster geometries proposed by Betts 
(see Fig. ^| to study the antiferromagnetic transition at 
half filling and generalize Betts' arguments to generate a 
set of clusters appropriate to study d-wave superconduc- 
tivity. 

To illustrate that the DCA recovers the correct re- 
sult as the cluster size increases, we plot in Fig. the 
DCA results for the Neel temperature Tn at half-filling 
as a function of the cluster size N c . Tn decreases slowly 
with increasing cluster size N c . As spin-correlations de- 
velop exponentially with decreasing temperature in 2D, 
the N c > 4 data falls logarithmically with N c , consis- 
tent with Tn = in the infinite size cluster limit. Thus, 
the Mermin- Wagner theorem is recovered for N c — > oo. 
The clusters with N c = 2 and N c — 4 are special be- 
cause their coordination number is reduced from four. 
For N c = 2 the coordination number is one and hence 
a local singlet is formed on the cluster for temperatures 
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FIG. 1: Cluster sizes and geometries used in our study. 
The shaded squares represent independent d-wave plaquettes 
within the clusters. In small clusters, the number of neigh- 
boring d-wave plaquettes Zd listed in table Q is smaller than 
four, i.e. than that of the infinite lattice. 
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FIG. 2: Neel temperature at half-filling when U — 8t versus 
the cluster size. Tn scales to zero in the infinite cluster size 
limit. The solid line represents a fit to the function A/(B + 
ln(iV c /2)) obtained from the scaling ansatz £(Tjv) = L c . For 
N c — 2 a local singlet and for N c = 4 the RVB state suppress 
antiferromagnetism. 

below J ~ t 2 /U. In the N c = 4 site cluster, the coordi- 
nation is two, so fluctuations of the order parameter are 
overestimated and the resonating valence bond state Q 
is stabilized. Hence, antiferromagnetism is suppressed in 
these cluster sizes and their corresponding Tn does not 
fall on the curve. 

We now turn to the main focus of this Letter, i.e. the 
search for a possible KT instability to the superconduct- 
ing state. To check that the DCA formalism is able to 
describe such a transition, we first tested the DCA-QMC 
code on the negative U, i.e. attractive Hubbard model 
which is known to exhibit a KT instability to an s-wave 
superconducting state [10|. We find that the DCA in- 
deed produces a finite temperature s-wave instability to 
the KT superconducting state. Due to the local nature 
of the s-wave order parameter, the DCA results converge 
rather quickly with cluster size. The DCA values for T c 
agree with those recently obtained in finite size QMC 
simulations ^(J- I n addition, we checked that our DCA- 
QMC code reproduces the results of other DMFT codes 
when N c = 1, and those of finite size QMC codes when 
the coupling to the self-consistent host is turned off. 

To identify a possible KT transition in the positive U 
Hubbard model we calculate the d x 2_„2-wave pair-field 
susceptibility P d for the clusters N c = AA, 8 A, 16A, 16B, 
18A, 20A, 24A and 26A. In contrast to the s-wave or- 
der parameter in the attractive model, the d-wave order 
parameter is non-local and involves four bonds or sites. 
Thus, large size and geometry effects have to be expected 
in small clusters. Similar to the cluster grading scheme 
Betts developed for magnetic order, we can classify the 
different clusters according to their quality for d-wave or- 
der. At low temperatures, local d-wave pairs will form, 
but phase fluctuations of the pair wave-function prevent 



the system from becoming superconducting. Since the 
DCA cluster has periodic boundary conditions, each four- 
site d-wave plaquette has four neighboring d-wave pla- 
quettes. However, as illustrated in Fig. ^ m small clus- 
ters, these are not necessarily independent and the effec- 
tive dimensionality may be reduced. 

Fig. n shows the arrangement of independent d-wave 
plaquettes in the clusters used in our study and their 
corresponding number Zd is listed in table [I] In the infi- 
nite system, Zd = 4. The N c = 4 cluster encloses exactly 
one d-wave plaquette (zd = 0). When a local d-wave 
pair forms on the cluster, the system becomes supercon- 
ducting, since no superconducting phase fluctuations are 
included. Thus, the N c = 4 result corresponds to the 
mean-field solution. In the 8A cluster, there is room for 
one more d-wave pair, thus the number of independent 
neighboring d-wave plaquettes Zd — 1. Since this same 
neighboring plaquette is adjacent to its partner on four 
sides, phase fluctuations are replicated and hence overes- 
timated as compared to the infinite system. The situa- 
tion is similar in the 16B cluster, where only two indepen- 
dent (and one next-nearest neighbor) d-wave plaquettes 
are found (z.d = 2). In contrast, Zd = 3 in the oblique 
16A cluster. We thus expect d-wave pairing correlations 
to be suppressed in the 16B cluster as compared to those 
in the 16A cluster. With the exception of the 18 A cluster, 
where neighboring d-wave plaquettes share one site and 
thus are not independent, the larger clusters 20 A, 24A, 
and 26A all have Zd = 4 and are thus expected to show 
the most accurate results. Hence, as the number of in- 
dependent neighboring d-wave plaquettes, Zd, is reduced 
from four, phase fluctuations are replicated due to peri- 
odic boundary conditions and thus overemphasized, sup- 
pressing pairing correlations and consequently T c . Note 
that the effects of finite size energy levels on the pair- 
ing correlations were pointed out in QMC simulations of 
Hubbard ladders [l9j ]. 



TABLE I: Number of independent neighboring d-wave pla- 
quettes Zd and the values of T^ T and Tl ln obtained from the 
Kosterlitz-Thouless and linear fits of the pair-field suscepti- 
bility in Fig. [3] respectively. 
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Tl in /t 
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(MF) 0.046 


0.056 


8A 
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-0.014 


-0.006 


18A 
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-0.043 


-0.022 


12A 
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0.011 


0.016 


16B 
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0.010 


0.015 


16A 
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0.021L0.008 


0.025L0.002 


20A 
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0.019 


0.022 


24A 
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0.016 


0.020 


26A 
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0.020 


0.023 



4 



Fig-EHshows the temperature dependence of the inverse 
d-wave pair-field susceptibility, 1/Pd, in the 10% doped 
system. Since a proper error propagation is severely ham- 
pered by storage requirements, we obtain the error-bars 
shown on the 16A results from a number of indepen- 
dent runs initialized with different random number seeds. 
Error-bars on larger cluster results are expected to be of 
the same order or larger. The results clearly substantiate 
the topological arguments made above. 

As noted before, the N c — 4 result is the mean-field 
result for d-wave order and hence yields the largest pair- 
ing correlations and the highest T c . As expected, wc 
find large finite size and geometry effects in small clus- 
ters. When Zd < 4, fluctuations are overestimated and 
the d-wave pairing correlations are suppressed. In the 8A 
cluster where Zd = 1 we do not find a phase transition at 
finite temperatures. Both the 12A and 16B cluster, for 
which Zd = 2, yield almost identical results. Pairing cor- 
relations are enhanced compared to the 8A cluster and 
the pair-field susceptibility Pd diverges at a finite temper- 
ature. As the cluster size is increased, Zd increases from 
3 in the 16A cluster to 4 in the larger clusters, the phase 
fluctuations become two-dimensional and as a result, the 
pairing correlations increase further (with exception of 
the 18A cluster). Within the error-bars (shown for 16A 
only), the results of these clusters fall on the same curve, 
a clear indication that the correlations which mediate 
pairing are short-ranged and do not extend beyond the 
cluster size. 

The low-temperature region can be fitted by the KT 
form P d = Aexp(2B/(T - T c ) - 5 ), yielding the KT esti- 
mates for the transition temperatures given in ta- 
ble[l] We also list the values T^ m obtained from a linear fit 
of the low temperature region, which is expected to yield 
more accurate results due to the mean-field behavior of 
the DC A close to T c . For all clusters with z d > 3 
we find a transition temperature T c « 0.023t ± 0.002t 
from the linear fits. We cannot preclude, however, the 
possibility of a very slow, logarithmic cluster size depen- 
dence of the form T C (N C ) = T c (oo) + B 2 / {C + ln(N c )/2) 2 
where T c (oo) is the exact transition temperature. In this 
case it is possible that an additional coupling between 
Hubbard planes could stabilize the transition at finite 
temperatures. 

In summary, we have presented DCA/QMC simula- 
tions of the 2D Hubbard model for clusters up to N c = 32 
sites. Consistent with the Mermin- Wagner theorem, the 
finite temperature antiferromagnetic transition found in 
the N c — 4 simulation is systematically suppressed with 
increasing cluster size. In small clusters, the results for 
the d-wave pairing correlations show a large dependence 
on the size and geometry of the clusters. For large enough 
clusters however, the results are independent of the clus- 
ter size and display a finite temperature instability to 
a d-wave superconducting phase at T c w 0.023i at 10% 
doping when U = At. 
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FIG. 3: Inverse d-wave pair-field susceptibility as a func- 
tion of temperature for different cluster sizes at 10% dop- 
ing. The continuous lines represents fits to the function 
P d = Aexp(2B/{T - T c ) - 5 ) for data with different values 
of Zd- Inset: Magnified view of the low-temperature region. 
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